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A technique inspired on quantum annealing is proposed in order to obtain the classical ground 
state of a spin-glass by tracking the full wavefunction of a given system within the subspace of matrix 
product states (MPS), using the density matrix renormalization group (DMRG). The technique is 
exemplified within the problem of obtention of the classical ground state of an Ising spin glass on 
ladder geometries. Its performance is evaluated and related to the entanglement entropy. 
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I. INTRODUCTION 

Global optimization is one of the most challenging nu- 
merical tasks. Simulated thermal annealing (STA) has 
been, for more than twenty years, one of the most popu- 
lar general purpose tools. Its strength relies on its ability 
to escape metastable minima via thermal fluctuations. 
Simulated quantum annealing (SQA)i is a more recent 
algorithm which takes profit from quantum fluctuations 
for the same purpose. Both methods have a clear phys- 
ical motivation: thermal annealing has been carried out 
by metallurgists for millenia. On the other hand, the 
first evidence of the superiority of quantum over thermal 
fluctuations in order to find the global minimum of a real 
physical system was only obtained in 1999, in the context 
of quantum spin glasses^. 

Numerical implementations of SQA have relied heavily 
on the quantum- classical analogy described by the path 
integral Monte-Carlo method (PIMC). Standard simu- 
lated thermal annealing is applied to an enlarged sys- 
tem, which consists of several replicas of the original 
system in interaction. Slightly different versions of the 
method have succeeded in the obtention of the ground 
state of classical spin glasses^, atomic clusters^, travel- 
ing salesman problem^, kinetically constrained systems^ 
and small protein-folding problems^, among others. SQA 
has proved to work less efficiently than STA in other 
problems, such as 3-SAT— and some benchmark ID 
potentialsiii. 

Despite the success of the PIMC-SQA, there are rea- 
sons to look for diff'erent implementationsii. First of all, 
PIMC simulations must be carried out at finite temper- 
ature, with lower temperatures requiring a larger num- 
ber or replicas and, therefore, higher computational cost. 
Also, some quantum systems suffer from the sign prob- 
lem, or from problems with the Trotter break-up. Sam- 
pling difficulties are always a risk for annealing schemes, 
and ensuring ergodicity is often a highly non-trivial task. 

Advancing in this line. Green's function Monte-Carlo 
technique was attempted by Stella and Santoro^^. Its 
main disadvantage is the necessity of good trial vari- 
ational wavefunctions. Real time evolution of the full 
wavefunction has been implemented by Suzuki and 
Okadap^, making use of both exact methods and the 



time-dependent density matrix renormalization group 
(DMRG) algorithroii. For systems where DMRG is ef- 
ficient, loss of adiabaticity in the form of Landau-Zener 
level crossings is the main problem of this approach. 

The proposal of this work is to perform the anneal- 
ing on the full wavefunction, thus overcoming most of 
the difficulties associated with PIMC. But, as opposed to 
the Suzuki-Okada approach, a real time evolution is not 
needed either. Instead, the ground state of the full sys- 
tem is computed exactly for a high intensity of the quan- 
tum fluctuations. As these fluctuations are decreased, 
the ground state is subsequently updated, by making 
the necessary small modifications of the previous ground 
state. The full wavefunction of the ground state is stored 
in the form of a matrix product state (MPS), and is com- 
puted variationally using the DMRG. 

As a benchmark case-study we have selected the ob- 
tention of the classical ground state of an Ising spin- 
glass with couphngs uniformly distributed in [—1,1]. The 
quantum fluctuations are provided by a transverse mag- 
netic field. We have chosen the ladders for the underlying 
topology, since it is the simplest case which presents gen- 
uine frustration while retaining the quasi-lD character 
which is required for the DMRG to attain its maximum 
efficiency. Albeit the problem is in class P, we will show 
it to be complex enough to be considered non-trivial. 

This article is organized as follows. Section |TT] intro- 
duces our model hamiltonian, summarizing its general 
features as applied to our featured topology. It is fol- 
lowed, in section HTll by a detailed description of the pro- 
posed quantum wavefunction annealing (QWA) method. 
In section HVl we analyze the results of the numerical ex- 
periments we have conducted. We conclude, in section 
IVl with a discussion of the advantages and limits of the 
method, and its possible extensions. 



II. RANDOM ISING MODEL IN A 
TRANSVERSE FIELD ON LADDERS 



Our model hamiltonian is the random Ising model in a 
transverse field (RITF), which we define on an arbitrary 
graph Q with N sitesi^: 
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where denotes nearest neighbours m Q, the J^.j are 
independent random variables uniformly distributed in 
[—1, 1], r is the transverse field and erf and erf are Pauli 
matrices at site i. We will work on the basis of eigen- 
states of cr^. Unless explicitly stated, we will assume the 
longitudinal fields /if to be zero. The first term in eq. [1] 
may be considered as a potential energy, and the second 
as a kinetic term: the magnitude of T measures the inten- 
sity of the quantum fluctuations. For F = 0, the ground 
state is composed of only two states, related by a trivial 
+ <^4> — symmetry, which we term the classical ground 
states (CGS) of the spin glass. As F — > oo, on the other 
hand, all spins should be pointing in the X-direction, 
therefore obtaining the following ground state: 
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All the wavefunction components of the ground state take 
the same value in our basis for F ^ oo. 

The behaviour of this system at T = is rather well 
known in I D^^i^^ thanks to an insightful RG analysis, in 
2E>i^ and 3E>i^ using quantum Monte-Carlo, and in ran- 
dom graphs with fixed connectivity's making use of the 
DMRG. In all these cases, the system presents a quantum 
spin-glass transition (QSGT) at a finite value of F = Fc, 
where the energy gap vanishes. Above that value, the 
system is said to be in a quantum paramagnetic regime, 
whilst below it behaves as a quantum spin-glass. 

The problem of finding the global CGS for equation 
[1] at F = is known to be in class P for D = 1 and 
for D = 2 in the absence of longitudinal magnetic field 
h^. In the 2D case with ft,^ 7^ 0, the 3D case^ or in 
random graphs with fixed connectivity^^, the problem is 
NP-complete. 

Rectangular ladders constitute the simpest topology in 
which the system described by equation [1] presents gen- 
uine frustration. Their size is characterized by two num- 
bers: L X w, where L is the length and w is their width, 
or number of legs. Figure [T] shows a specimen with size 
5x2. The classical spin-glass with ± J couplings on these 
ladders has been analyzed by a number of authors^^i^l. 
We will now discuss the complexity of the energy land- 
scape for F = and a few basic characteristics of the 
nature of the QSGT. 

Minimization of the classical hamiltonian (F = 0) can 
be done using STA or PIMC-SQA. For the system sizes 
under consideration, both methods yield robust estimates 
for the CGS energy following the schemes described by 
Santoro and coworkers^. In order to analyse the com- 
plexity of the energy landscape, we relax the parameters 



f----f- 

I I 
I I 



1:1 



FIG. 1: A rectangular ladder of dimension 5x2. Coupling 
constants Jij are associated to links of the graph. In the exam- 
ple of the figure, dashed (continuous) lines represent negative- 
AFM (positive-FM) links, and the system is frustrated. 



FIG. 2; Some local minima of the classical hamiltonian for 
three 40 x 2 samples, obtained with a non-robust STA in 
order to sample the low energy configurations. 



until the system becomes non-robust, i.e.: it yields dif- 
ferent results in different runs. At this point, the anneal- 
ing processes provide us a series of sample configurations 
which constitute local minima of the energy. With this 
purpose, we have applied a STA algorithm with a mul- 
tiplicative scheme for f3, i.e.: P rf3, from Pq = 0.1 
to Pmax = 10^, with r = 1 -I- 10^'"' and 10"* steps per 
temperature. The results for three samples of a 40 x 2 
ladder are shown in figure [H Each column contains the 
energy values obtained after the procedure was repeated 
20 times, providing several different local energy minima 
(about 10). In all the cases shown, the lowest energy 
corresponds to the CGS. This result points to a complex 
energy landscape for the spin-glass ladders. 

The energy per spin of the CGS seems to converge, 
in the case of the 2-legged ladders, to e^^^ « —0.64. In 
absence of frustration, this value would be —0.75. In 
the case of the ladder with 4 legs, it converges to e''*-' ~ 
—0.71, with the unfrustrated value being —0.875. 

At T = 0, the system described by eq. [1] presents a 
quantum spin-glass transition (QSGT) at a finite value 
of F = Fc. The finite-size DMRG algorithm, suitably 
adapted for our case^S, can be used to characterize this 
transition. The behavior of some relevant observables 
has been traced in figure [3l (a) the energy gap AS, (b) 
the maximum block entropy Smax^ given by 



Smax = -TrpLln(/9L) 



(3) 



where pL is the reduced density matrix for the left half 
of the ladder, and (c) the spin glass susceptibility xsGj 
as defined by the following formula: 
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FIG. 3: Behavior of the energy gap A_B, the maximum block 
entropy Smax and the spin-glass susceptibility XSG for a sam- 
ple spin-glass ladder of dimension 40 x 2, as a function of F. 
The QSGT is marked very clearly by both xsG and AE. On 
the other hand, Smax seems to have a more erratic behavior. 



XSG = 




i.e.: a small longitudinal magnetic field applied at 
site j, generates a magnetization response on each site 
i, which is measured (and squared, so as to disregard its 
sign) ; the results are summed over all sites i and averaged 
over all sites j. The value of Fc changes from sample to 
sample. Within a single sample, AE vanishes and xsG 
diverges at the same value of Fc, as it is shown in figure 
[31 For that sample. Fc = 0.6. The entropy is seen to 
present a more complex behavior. 

Some insight can be gained by tracing a few individ- 
ual components of the full wavefunction, which is pos- 
sible within the DMRG framework^^. Fi gure |4] shows 
that, well within the paramagnetic regime, all wavefunc- 
tion components take the same value. As we decrease F, 
the configurations with low energy start to increase its 
weight in the ground state wavefunction, while the con- 
figurations with high energy start to decrease. At the 
critical point, all the configurations have started their 
decay, except the one with the lowest energy. 

More theoretical and analytical work is needed in order 
to fully characterize this QSGT. Some comments about 
it are provided in section [V] This preliminary analysis 
can be summarized as: (a) the energy landscape of the 
classical problem is complex, and (b) there is a quantum 
phase transition at a finite value of F. 
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FIG. 4: Wavefunction components for some selected config- 
urations, as a function of F, for a sample ladder with size 
40 X 2. All the configurations but one are local minima of the 
classical energy function, obtained using a non-robust STA. 
The transition point is marked with an vertical line. 



III. QUANTUM WAVEFUNCTION 
ANNEALING OF SPIN GLASS SYSTEMS 

Quantum wavefunction annealing (QWA) is an imple- 
mentation of simulated quantum annealing (SQA) which 
does not require real time evolution, configurations sam- 
pling or Trotter splitting, ft makes extensive use of ma- 
trix product states (MPS) ^^i^^ and the density matrix 
renormalization group (DMRG)^. For the benefit of the 
readers which are not acquainted with those techniques, 
the first paragraph of this section describes a simpler for- 
mulation based on the Lanczos algorithm, but which can 
only be used for very small lattices. 

The Lanczos diagonalization method improves its per- 
formance dramatically if a good seed to the real ground 
state of the system is provided^^. The hamiltonian given 
by equation [T] has a simple behaviour for F ^ oo: its 
ground state is known to be given by the state \X) de- 
fined in equation[5) Therefore, starting with such a state 
as a seed, it is fairly easy for the Lanczos algorithm to 
obtain the exact ground state for a finite but high value 
Fq. This value may be decreased at small finite steps, 
F ^ F — AF, and the ground state for the previous value 
of F can be employed as the seed for the computation of 
the ground state at the new value, fn quantum annealing, 
the adiabatic theorem is needed in order to prove that 
convergence to the CGS is sure in the limit AF 0. In 
our case, a weaker result is enough: convergence is sure 
as long as there is a finite overlap between the ground 
states at any two consecutive values of F. 

The obvious drawback of the previous algorithm is 
the size of the Hilbert space, dim(7i) = 2^. Using 
the plain Lanczos algorithm, all the wavefunction com- 
ponents should be stored and acted upon. A solution 
is to choose a low-dimensional subspaec which is known 
to contain the ground state of the hamiltonian given by 
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equation [T] for all values of T. Matrix product states 
(MPS) can provide such a subspace in some cases. These 
states may be written down as: 

|vp) =^tr(A^iA^^...A^")|si,S2,... (5) 

i 

where |si, S2, • • • , sn) are the eigenstates of and the 
A'^' are 2N matrices of dimension m x m, whose entries 
are the variational parameters of our Ansatz. The di- 
mension of the MPS subspace is, therefore, bounded by 
2Nm?. Both the computational cost and the accuracy of 
the method depend strongly on the dimension m of the 
matrices, whose physical meaning is the following: for 
all possible left-right splittings of the system, the ground 
state is approximated retaining m states to represent the 
left part and other m states for the right part. For any 
possible state of the system, its representation as an MPS 
becomes exact for m large enough. Of course, the dimen- 
sions m can be made local, m^, if necessary. 

The DMRG is a variational method within the sub- 
space of the MPS2^. It profits from the use of a density 
matrix in order to select the m states for the left and 
right blocks which fit best to our current approximation 
of the global ground state. The neglected eigenvalues of 
the density matrix provide a way to monitor the accuracy 
of the method. In all our applications, the tolerance to 
the total neglected probabilities is fixed beforehand to a 
certain value 77, and the number of retained states m{ri), 
is adapted in consonance. 

The MPS represent faithfully ground states of local ID 
hamiltonians^. In other cases, at this stage, only edu- 
cated experience can decide whether they are appropri- 
ate or not. They are specially suited to ID and quasi-lD 
problems (e.g. ladders, trees), although natural exten- 
sions to higher dimensions are in active devclopmen t^^'^^ . 
A crucial tool is von Neumann's block entropy, S = 
— tr(plog2 /o). In crude terms, m should scale as exp(5) 
in order to obtain an accurate DMRG method. In ID, 
the entropy S is known to be bounded for non-critical 
systems and to scale as log(A^) for a critical one^. For 
higher dimensions, the area law predicts the entropy to 
scale as L'^~^ out of criticality, where L is a typical di- 
mension of the systen*^. 

As in the Lanczos case, the DMRG can benefit 
from a good seed when searching the ground state 
of a hamiltonian through the use of the wavefunction 
transformation^, which allow to use the solution of an 
RG-step as an initial step for the next one. Having found 
the ground state of a certain hamiltonian H , a few finite- 
size sweeps will adapt it to become the ground state of a 
slightly modified hamiltonian H' . In normal cases, this 
procedure is not needed, since the number of sweeps re- 
quired for convergence to the ground state for any hamil- 
tonian is small. It may become very useful when (a) it is 
known that the ground state of a certain hamiltonian can 
be represented as a MPS with low m but (b) the proba- 
bility of the DMRG getting stuck at an excited state is 



very high. This is the case for the ground state of[l]in 
the quantum spin-glass phased. 

Our proposed QWA algorithm is, therefore, the follow- 
ing one: 

(a) The ground state of the systcm[T]is obtained for very 

high r. 

(b) The transverse field V is decreased F ^ F — AF. 

(c) A few finite-size sweeps of the DMRG are carried 

out, which adapt the ground state to the current 
value of F. 

(d) Go to (b) if F is not yet zero. 

(e) Measure the energy and classical ground state. 

This approach is deterministic, i.e.: not limited by 
sampling problems. It works directly with the quan- 
tum hamiltonian, therefore does not require any Trot- 
ter break-up, and works directly at T = 0. Since it 
does not simulate real time evolution, it is not prone to 
Landau-Zener level-crossings. Loss of adiabaticity is not, 
therefore, a crucial issue. As it was stated, it is enough 
to ensure that the overlap between the ground states at 
consecutive values of F is finite. A high value of AF is al- 
lowed far away from the quantum phase transition, while 
a more reduced value will be taken near it. Our precise 
adaptive reduction schedule is described in the following 
section. 

On the negative side, it is based on a method which 
is specially designed for ID systems. We will say more 
about this in the conclusions. Its main weakness stems 
from the inability of the MPS to represent faithfully the 
ground state for all values of F for higher dimension. The 
efficiency and practical issues of the implementation are 
analyzed in the next section. 

IV. NUMERICAL EXPERIMENTS 

In order to study the applicability of QWA, we have 
generated random samples of ladders with w = 2 and 
w = A legs, and lengths ranging from L = 20 to L = 160. 
We have also considered random graphs of fixed connec- 
tivity K = Z and, as a check, we have also analysed 
linear chains (although the optimization problem is triv- 
ial in that case). The annealing scheme has always been 
the same: starting with Fq = 3 and using an adaptative 
reduction, F ^ F - AF with 

AF = min(0.5,0.1/S') (6) 

where S is the maximum block entropy for the previous 
value of F. This way we ensure that, near the QSGT - 
when the entropy is higher- the steps are shorter. We in- 
sist in this fact: the adiabatic theorem and Landau-Zener 
level crossings are not the limitations of this method. 



5 



Therefore, a decrease in the size of the steps does not 
necessarily increase the probabihty of success. 

The last annealing step was always taken with Tmin = 
0.01, a value which is low enough for practical purposes. 
The obtcntion of the minimum energy configuration was 
carried out measuring the z-component of the spin po- 
larization of all sites after convergence for F = Tmin- Of 
course, due to the obvious symmetry + —, the expec- 
tation value for each of these components is zero. We 
have used a common procedure in quantum spin-glass 
calculations: the insertion of a very small longitudinal 
magnetic field /i^ = 10"^ in a single site, selected at ran- 
dom, in order to break the symmetry. In the QSG phase, 
because of the divergence of the spin-glass susceptibility 
discussed above, an infinitesimal localized magnetic field 
polarizes the full sample. 

Table U provides the basic set of results. For each ge- 
ometry we have selected 20 samples and performed the 
QWA algorithm on them. We have also chosen differ- 
ent values of the neglected probabilities tolerance, rj. For 
each run, QWA is said to obtain a success if its mini- 
mum energy is equal to the value obtained by STA and 
PIMC-SQA. QWA has had success, in all the attempted 
geometries, whenever rj = 10~^. When the tolerance was 
decreased, the CGS was missed with higher probability, 
although the method can be seen to be rather robust for 
spin-glass ladders, since the tolerance must be raised up 
to lO""^ in order to decrease the probability of success to 
50%. 

We have also tried to check the method with a different 
graph topology: random graphs with fixed connectivity, 
K = 3. For = 20 sites, employing a tolerance of 
rj — 10^^, the method provides again 100% of success. If 
the size is increased, that tolerance can not be set, since 
the number of states m to be kept becomes prohibitive. 
Therefore, we have carried out the experiments with rj = 
10"^ for aN = 100 sample, and the probability of success 
reduces to 45%. 

The time for the QWA algorithm scales as a power law 
of the system size: T w L". For linear chains, using the 
results from table |ll a w 1.3 ± 0.1. For 2-leggcd ladders, 
a « 2±0.1. In the case of the 4-legged ladders, the power 
law fit has a higher errorbar, with a surprising exponent 
of asa 1.5 ±0.2. 



Geom. 




Success 


Time 


20 


10"* 


100% 


1.2 ± 0.2 


40 


10"* 


100% 


3.9 ±0.8 


80 


10"* 


100% 


10 ±2 


160 


10"* 


100% 


21 ±4 


320 


10"* 


100% 


56± 10 


20 X 2 


10"* 


100% 


17 ±8 


40 X 2 


10"* 


100% 


60± 15 


80 X 2 


10"* 


100% 


240 ± 30 


160 X 2 


10"* 


100% 


600 ± 110 


40 X 2 


10"^ 


80% 


26 ±3 


40 X 2 


10"3 


50% 


13 ± 1 


20 X 4 


10-* 


100% 


1600 ± 600 


40 X 4 


10"* 


100% 


6800 ± 2000 


80 X 4 


10"* 


100% 


14000 ± 3000 


160 X 4 


10"* 


100% 


27000 ± 3000 


RG-20 


10-* 


100% 


190 ± 150 


RG-100 


lo-'* 


45% 


15000 ± 7000 



TABLE I: Numerical results for the QWA method. The first 
column states the geometry of the system under study: L 
if it is a linear chain, L x w if it is a ladder, and RG-A^ if 
it is a random graph with connectivity K — 3. The second 
gives the tolerance for neglected probabilities in the DMRG, 
rj. The third provides the percentage of success. The fourth 
column provides the average time (in seconds) for the QWA 
of a single sample. 

and therefore we may expect T w Lexp(25c) ~ 
Thus, the theoretical prediction is a ~ 2a ± 1. 

In the linear chain case, the entropy at criticality was 
predicted by Rcfacl and Moore^ to have a coefficient 
a = ln(2)/6 w 0.11, in agreement with our own numeri- 
cal measurements. Therefore, the theoretical prediction 
for the a exponent is « 1.22, which is not far from the 
1.3 ±0.1 obtained numerically. In the case of the 2-legged 
ladder, there is no theoretical estimate, but our own nu- 
merical simulations provide a value a ~ 0.55 ±0.1, thus 
giving an estimate for a around 2.1, in agreement with 
the numerically observed value a « 2 ± 0.1. 



A. QWA and block entropy 

This polynomial growth can be theoretically explained. 
For quasi-lD systems at criticality, the von Neumann en- 
tropy grows as a logarithm of the system size^, Sc{L) w 
alog(L)±6. The value of a has been related, in some ran- 
dom systems, to the central charge of the associated con- 
formal field theory^. On the other hand, the maximum 
number of retained states scales as the exponential of the 
entropy, m exp(S'), and the time for a DMRG sweep 
scales as T « Lrri^. The most expensive DMRG sweep 
for a QWA simulation takes place at the critical point. 



V. CONCLUSIONS AND FURTHER WORK 

In this paper we have introduced a method for the 
obtention of the minimum energy configuration of a 
spin-glass, based on the annealing of the full wavefunc- 
tion represented as a matrix product state (MPS) using 
the density matrix renormalization group (DMRG). The 
method has been termed quantum wavefunction anneal- 
ing (QWA), and we have assessed and quantified its ef- 
ficiency, by comparing its results with those provided by 
other robust methods. For spin-glass ladders of 2 and 4 
legs with lengths ranging from 20 to 160, the method has 
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always provided the optimum solution when the proba- 
bility tolerance is set to ?/ = 10^*. The running time 
scales as for 2-legged ladders, and, surprisingly, as 
Jqj. 4_iegged ones. This anomaly requires further 
clarification, perhaps in the line of thought of Ferraro et 
al^. For random graphs with fixed connectivity K = 3, 
on the other hand, the tolerance has to be increased in 
practice to 10~^, and the probability of success falls to 
45% for N = 100. 

The most crucial parameter which determines the suc- 
cess of the method is the probability tolerance rj. When- 
ever we were able to set rj = 10~®, for whichever topology, 
the system always attained the optimum configuration. 
Unfortunately, for topologies which are not quasi-lD, the 
number of retained states in DMRG, to, grows very fast 
when the tolerance is decreased. Other parameters, such 
as the annealing velocity, do not have such a direct rel- 
evance to the quality of the results. Loss of adiabaticity 
is not the main issue for this method: a finite overlap 
between the ground state wavefunctions at consecutive 
annealing steps is enough to ensure convergence, as long 
as the number of retained states to is large enough for our 
Ansatz wavcfunction to represent both of them faithfully. 

Therefore, it seems reasonable to think that QWA can 
be converted into a method to find the optimum con- 
figuration of all quasi-lD systems in polynomial time. 
This suggests that, in fact, no quasi-lD problem can be 
NP-complete. On the other hand, for problems which 
are known to be NP-complete, such as the random Ising 
spin-glass in a 2D lattice with a longitudinal field, or on 
random graphs of fixed connectivity, the time for an algo- 



rithm which is able to obtain the optimum with certainty 
should scale exponentially with the size of the system. 

This fact puts a limit to possible extensions of the 
QWA algorithm described in this paper. The natu- 
ral extension of the MPS are the tensor product states 
(TPS) analyzed by Nishino and coworkers^ or the pair- 
entangled particle states (PEPS)^ of Verstraete and 
Cirac. The computational power of this type of states 
has been analyzed recently^^, pointing to the fact that 
any attempt to solve NP-complete problems with an ex- 
act QWA algorithm using these states would require an 
exponential time, which is compatible with the usually 
believed notion that P 7^ NP. 

The most promising lines for future work are, therefore, 
in the development of heuristic algorithms, running in 
polynomial time, which may give the absolute minimum 
energy state with high probability for some problems. A 
possibility is the development of a QWA algorithm using 
TPS (or PEPS) with a fixed number of retained states 
TO. Another possibility is to use the DMRG approach to 
non-equilibrium classical problems^ to develop a classi- 
cal wavefunction annealing algorithm, which would solve 
the Fokker-Planck equation associated to simulated ther- 
mal annealing. 
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